Introducción

Si es su primera vez trabajando con los datos, no está muy familiarizado con el paquete o simplemente quiere conocer más sobre el proyecto y el objetivo de estos ejemplos y el paquete verdata, consulte: https://github.com/HRDAG/CO-examples/blob/main/Introducción/output/Introducción.html antes de continuar.

En este ejemplo, se ilustrará el proceso para obtener los datos observados imputados y la estimación por año del subregistro de víctimas de desaparición forzada. Específicamente, se replicará la tabla 9 anexo metodológico del proyecto.

Autenticando e importando la base de datos (réplicas)

Se comienza autenticando e importando la base de datos de desaparición, esto a través de dos funciones del paquete verdata: las funciones confirm_files y read_replicates. La autenticación de los datos es pertinente dado que estos fueron publicados con la licencia de atribución 4.0 internacional de Creative Commons (CC BY 4.0). Esta licencia permite la distribución y modificación de la información. Considerando que usted pudo haber llegado a estos datos por medio de diferentes fuentes, es importante que sepa si han sido modificados o no, para lo que puede hacer uso de estas dos funciones.

La función confirm_files autentica los archivos que han sido descargados. Considerandoque cada violación tiene hasta 100 réplicas, esta función permite autenticar cada uno de estos archivos sin necesidad de leerlos a R. Esto, en caso de querer ahorrar recursos computacionales, o en caso de que no vaya a realizar su análisis con todas las réplicas. Esta función devolverá una tabla con dos columnas: una indicando la ruta del archivo y otra indicando si el archivo es igual al publicado. En caso de que al menos uno de los archivos no sea igual, la función devuelve el mensaje “Some replicate file contents do not match the published versions”.

confirmar <- verdata::confirm_files(here::here("verdata-parquet/reclutamiento"),
                                    "reclutamiento", c(1:10))

Además, la función read_replicates permite 2 cosas: leer las réplicas a R en una sola tabla (ya sea a partir de un formato csv o parquet) y verificar que el contenido de las réplicas sea exactamente igual al publicado. Cuando el argumento crash tiene su valor por default (TRUE), la función retorna un objeto (data frame) si el contenido es igual, y el mensaje “The content of the files is not identical to the ones published. This means the results of the analysis may potentially be inconsistent.” si el contenido de la base fue previamente alterado/modificada, lo que quiere decir que los análisis que el usuario realice hacer serán inconsistentes y llevarán a resultados erróneos. Este último error significa que nos datos no se han leído a R. Si por alguna razón, usted quiere leer los datos a pesar de saber que no son los mismos datos originamente publicados, puede cambiar el argumento crash a FALSE, y, en ese caso, podrá ver los datos, junto con el mismo mensaje de advertencia.

replicas_datos <- verdata::read_replicates(here::here("verdata-parquet/reclutamiento"),
                                           "reclutamiento", c(1:10))

paged_table(replicas_datos, options = list(rows.print = 10, cols.print = 5))

Vemos que tenemos 192 260 registros, nuestras réplicas van desde la número 1 hasta la 10. Además, nuestros datos tienen información sobre la categoría de edad de la víctima, el presunto perpetrador, el sexo, el año del hecho, la pertenencia étnica, entre otros. Sin embargo, para centrarnos en un análisis más específico, tal como el realizado para el informe metodológico, procederemos a transformar y/o filtrar algunas variables.

Filtrando las réplicas acorde con el filtro del informe metodológico

La función filter_standard_cev nos permite transformar o filtrar la información. Por ejemplo, para el caso de reclutamiento optamos por analizar a las víctimas menores de 18 años, es decir, a aquellas víctimas que pertenecen a la categoría de infancia y adolescencia, generando la categoría “MENOR”. Adicionalmente aquellos menores de edad que se documentaron como víctimas de la ex-guerrilla FARC-EP en años posteriores a 2016 pasaron a ser víctimas de otras guerrillas, ya que este primer grupo oficialmente dejó de existir después de dicho año (perp_change = TRUE).

replicas_filtradas <- verdata::filter_standard_cev(replicas_datos,
                                                   "reclutamiento", 
                                                   perp_change = TRUE)

paged_table(replicas_filtradas, options = list(rows.print = 10, cols.print = 5))

Víctimas documentadas

Después de aplicada la anterior función es momento de obtener la información de las víctimas documentadas por sexo. Estos datos en particular son aquellos que ya venian en la base integrada y que en ocasiones contenían campos faltantes en algunas de las variables. Usaremos la función summary_observed para calcular dicha documentación.

Como se puede ver, los argumentos de la función son: 1) la violación a analizar reclutamiento; 2) los datos replicas_filtradas; 3) strata_vars, que para este ejemplo será la variable de yy_hecho, porque estamos analizando la desaparición por año; 4) le sigue el argumento de conflict_filter que filtra a aquellas personas que fueron víctimas dentro del marco del conflicto armado (variable is_conflict == TRUE) o no (variable is_conflict == FALSE). Para el hecho de reclutamiento todas las víctimas fueron reclutadas dentro del conflicto armado.

Esta función también incluye un argumento denominado 5) forced_dis_filter, que aplica únicamente a esta violación. Esta indica si la víctima fue desaparecida de forma forzada (forced_dis == TRUE) o no (forced_dis == FALSE). Para otras violaciones este argumento siempre será FALSE.

También contamos con otros argumentos: 6) edad_minors_filter que filtra por víctimas menores de edad (edad_minors_filter = TRUE) documentadas por los proyectos y/o instituciones; 7) include_props que permite incluir el cálculo de las proporciones para las variables de interés (include_props = TRUE).

tabla_documentada <- verdata::summary_observed("reclutamiento",
                                               replicas_filtradas, 
                                               strata_vars = "sexo",
                                               conflict_filter = TRUE,
                                               forced_dis_filter = FALSE,
                                               edad_minors_filter = TRUE,
                                               include_props = TRUE)

paged_table(tabla_documentada, options = list(rows.print = 4, cols.print = 4))

Posterior a esto se procede a aplicar la función de combine_replicates que, como su nombre lo indica, permite combinar las réplicas para obtener lo que denominamos “la media de la estimación puntual” junto con el intervalo de confianza que permite dimensionar la incertidumbre de la imputación. Para esta función se siguieron las reglas de combinación de Rubin, que, si desea estudiar con más detalle de qué se trata, el libro Flexible Imputation of Missing Data de Stef van Buuren aborda paso a paso este proceso.

Ahora, esta función se compone de los siguientes argumentos: la violación a analizar desaparicion; 2) tabla_documentada, es decir, el data frame derivado de la función summary_observed; 3) la base de datos filtrada replicas_filtradas; 4) strata_vars que será nuevamente la variable de año del hecho; 5) conflict_filter que filtra a aquellas personas que fueron víctimas dentro del marco del conflicto armado (variable is_conflict == TRUE) o no (variable is_conflict == FALSE).

Esta función también incluye un argumento denominado 5) forced_dis_filter que cual aplica únicamente a la violación de desaparición. Esta indica si la víctima fue desaparecida de forma “forzada” (forced_dis == TRUE) o no (forced_dis == FALSE). Para otras violaciones este argumento siempre será “FALSE”.

También contamos con otros argumentos: 6) edad_minors_filter que filtra por víctimas menores de edad (edad_minors_filter == TRUE); 7) include_props que permite incluir el cálculo de las proporciones para nuestras variables de interés (include_props == TRUE); y 8) digits que es un argumento opcional con el cual se puede establecer el número de dígitos para redondear los resultados (que por defecto es 2).

tabla_combinada <- verdata::combine_replicates("reclutamiento",
                                                tabla_documentada,
                                                replicas_filtradas, 
                                                strata_vars = "sexo",
                                                conflict_filter = TRUE,
                                                forced_dis_filter = FALSE,
                                                edad_minors_filter = TRUE,
                                                include_props = TRUE)

tabla_combinada <- tabla_combinada %>% 
  select(sexo, observed, obs_prop_na, obs_prop, imp_lo, imp_mean, imp_hi,
         imp_lo_p, imp_mean_p, imp_hi_p )

paged_table(tabla_combinada, options = list(rows.print = 10, cols.print = 5))

La primera columna se refiere a las víctimas que fueron registradas en las bases de datos como menores de edad, y que, sin embargo, no conocemos el sexo de 2 355 víctimas; y como se indicó, en esta columna excluímos víctimas mayores de 18 años y personas de las que no conocemos la edad. Estas últimas fueron objeto de la imputación estadística múltiple y por tal razón vemos que en la columna imp_mean hay una mayor cantidad de víctimas: primero porque esos 2 355 se imputaron y por ende estas víctimas “pasaron” a algunas de las categorías de sexo; y segundo, porque también incluímos a esas víctimas de las que no conocíamos la edad. Es decir, estas nuevas columnas incluyen todas las víctimas, independientemente si pasaron por un proceso de imputación o no.

Así, después del proceso de imputación estadística y con un con un nivel de confianza del 95% se evidencia que hubo entre 10 977 y 11 651 víctimas hombres, con un promedio de 11 314. Es decir, esto implica que este promedio es la mejor estimación puntual respecto al número de víctimas de esta categoría, no obstante, hay que tener en cuenta que siempre tendremos la incertidumbre de la imputación y que dicho fenómeno estará representado por el intervalo de confianza del 95%. Procedemos a estratificar para el posterior proceso de estimación del subregistro de víctimas.

Proceso estratificación para estimaciones

Con el fin de controlar la heterogeneidad en las probabilidades de captura (ver más de este concepto en el anexo metodológico del proyecto) se estratifica la información de acuerdo al análisis a realizar. En este caso, como queremos estimar el subregistro de reclutamiento por sexo, se estratifica por etnia, sexo y periodo presidencial. Es importante que usted como usuario vea que este proceso es netamente artesanal, es decir, usted puede usar su propio código o funciones para realizar este proceso que, en nuestro caso, será a través de una función previamente creada (fuera del paquete verdata) para facilitar este ejercicio:

stratify <- function(replicate_data, schema) {
    
    schema_list <- unlist(str_split(schema, pattern = ","))
    
    grouped_data <- replicate_data %>%
        group_by(!!!syms(schema_list))
    
    stratification_vars <- grouped_data %>%
        group_keys() %>%
        group_by_all() %>%
        group_split()
    
    split_data <- grouped_data %>%
        group_split(.keep = FALSE)
    
    return(list(strata_data = split_data,
                stratification_vars = stratification_vars))

}

Entonces, en primera instancia creamos una función que necesita de dos argumentos:

• El argumento replicate_data se refiere a un data frame a estratificar, que en nuestro caso es replicas_estratos, es decir, la información con las nuevas variables periodo_pres y etnia2.

• El segundo argumento son las variables de estratificación (schema). Recordemos que la estratificación es un instrumento para aislar la heterogeneidad, entonces estas son variables que pensamos pueden afectar la probabilidad de registro de las víctimas y, por lo tanto, queremos agrupar las víctimas con características similares.Todas estas variables deben encontrarse en el objeto replicas_estratos.

En términos generales, lo que hace esta función es: primero agrupa por las variables de estratificación y guarda en una lista llamada strata_data esta información. En ese sentido, cada elemento de la lista es una tabla con las víctimas que hacen parte de ese estrato. En segundo lugar, se define el nombre de cada estrato, para poder identificarlos cuando estimemos. Para esto se retorna una lista llamada stratification_vars que contiene las combinaciones de las variables, es decir, el nombre del estrato.

A continuación se aplica la función:

schema <- ("replica,periodo_pres,etnia2,sexo")

listas <- stratify(replicas_filtradas, schema)

replica6_grupo6 <- listas[["strata_data"]][[6]]

paged_table(replica6_grupo6, options = list(rows.print = 10, cols.print = 5))

replica6_grupo6_estrato <- listas[["stratification_vars"]][[6]]

paged_table(replica6_grupo6_estrato, options = list(rows.print = 10, cols.print = 5))

Estimación víctimas por año del hecho

Ahora, definidos los estratos, se calculan las estimaciones para esta violación en particular; para este paso se usa la función mse del paquete verdata. Esta función permite preparar los datos, revisar si hay estimaciones precalculadas para ese estrato y estimar, en caso de que no las haya. Como se indicó al principio, esta función tomará como insumo la información de las fuentes, es decir, aquellas columnas que comienzan por in_. También filtra por fuentes válidas, es decir, las fuentes que cuentan con al menos una víctima en ese estrato. Para que un estrato sea estimable, se requiere un mínimo de 3 fuentes válidas. Si un estrato no es estimable la función arrojará NA.

Como se mencionó, considerando que el proceso de estimación toma tiempo y recursos computacionales, esta función le permite usar estimaciones ya calculadas en el proyecto. Para esto, usted debe descargar las estimaciones publicadas acá.

mse(
  stratum_data,
  stratum_name,
  estimates_dir = NULL,
  min_n = 1,
  K = NULL,
  buffer_size = 10000,
  sampler_thinning = 1000,
  seed = 19481210,
  burnin = 10000,
  n_samples = 10000,
  posterior_thinning = 500
)

Algunos de los argumentos de esta función se explican de la siguiente forma¹

• stratum_data: Data frame que incluye la información del estrato de interés (data frames que creamos antes).

• stratum_name: Es el nombre del estrato.

• estimates_dir: Es la ruta (opcional) o el path de la carpeta o archivo que contiene las estimaciones precalculadas. Esto le permite a la función buscar entre las estimacines si el estrato que usted quiere analizar ya fue estimado.

estimaciones_dir <- here::here("estimaciones")

Al final el resultado será un data frame con cinco columnas: una columna que indica si el estrato es válido o no (TRUE o FALSE); el número de muestras N de la distribución posterior (NA si el estrato no es válido); las fuentes válidas que se usaron en la estimación valid_sources; el número de observaciones de las listas que son válidas del estrato (n_obs) y el nombre del estrato stratum_name. Si un estrato es estimable, este dataframe va a contener 1000 muestras para cada estrato. Esto lo ilustramos a continuación.

estimacion <- purrr::map2_dfr(.x = listas$strata_data,
                              .y = listas$stratification_vars,
                              .f = verdata::mse,
                              estimates_dir = estimaciones_dir)

# Separar variables
estimacion <- estimacion %>% 
    mutate(stratum_name = paste(pull(stratum_name, 1),
                                pull(stratum_name, 2),
                                pull(stratum_name,3),
                                pull(stratum_name,4),
                                sep = "-")) %>% 
    separate(stratum_name, into = c("replica", "periodo_pres", "etnia2", "sexo"), 
             sep = "-", remove = FALSE, extra = "merge") %>% 
    rename(replicate = replica) %>% 
    select(validated, N, valid_sources, n_obs, stratum_name, replicate, periodo_pres, etnia2, sexo)
    


paged_table(estimacion, options = list(rows.print = 10, cols.print = 8))

Estos son los resultados de mse. Vemos las cinco (5) primeras columnas mencionadas en el que cabe destacar que todos los estratos son válidos, es decir, no vemos ningún NA en nuestras columnas; lo que indica que en cada estrato hay al menos 3 listas con al menos 1 víctima. Además, vemos nuestro N el cual evidencia 1000 muestras aleatorias de la distribución posterior de la cantidad de víctimas probables para cada estrato y réplica (por ejemplo, 1000 números (o resultados) en N para el estrato 1990_1993-ETNICO-HOMBRE para la réplica 1):

union_mse_summary <- estimacion %>% 
    filter(replicate == "R1") %>% 
    group_by(stratum_name) %>% 
    summarize(N_estimaciones = n())

paged_table(union_mse_summary, options = list(rows.print = 10, cols.print = 5))

Análisis de la distribución posterior

Ahora, lo que sigue - que es un paso opcional - es análizar estos resultados desde el punto de vista visual, esto a través de una gráfica de densidad la cual nos va a permitir ver - gráficamente - la muestra de la distribución posterior para N. Para esto tomaremos -por ejemplo- la réplica 1 (R1) y el estrato 1998_2001-ETNICO-HOMBRE:

tabla_grafica <- estimacion %>% 
    filter(replicate == "R1") %>% 
    filter(stratum_name == "R1-1998_2001-ETNICO-HOMBRE") 

densidad_graph <- tabla_grafica %>% 
    ggplot() +
    geom_density(aes(x = N), color = "black") +
    theme_minimal() 

densidad_graph

Vemos varias cosas: lo primero es que tuvimos muestras menores a 600 y algunas de las muestras fueron mayores a 1000. La incertidumbre es menor hacia la parte izquierda de la distribución, ya que el número minimo que estimamos con LCMCR es el número documentado. Es decir, no es posible ni lógico que los resultados evidencien estimaciones menores al número documentado, que para este estrato en particular es 423. Esta situación es opuesta para el límite superior el cual puede ser mucho mayor (pero no mayor a la población colombiana). Por esta razón en este y otros ejemplos veremos que el límite superior (que lo abordaremos en la próxima sección) presenta más variación (y por tanto incertidumbre) que el límite inferior, es decir, gráficamente y analíticamente hablando nuestras distribuciones no siempre serán simétricas.

Otra característica que podemos observar es que no hay multimodalidad en este estrato, es decir, no vemos dos o más picos en nuestra distribución (es unimodal), por lo que podemos inferir que este es un estrato con heterogeneidad controlada. No obstante, vemos que hay varianza (pero no considerable), es decir una mayor -pero no exagerada- dispersión de las muestras.

Combinación de las estimaciones

Nuestro último paso es combinar nuestras estimaciones. Esto lo haremos a través de nuestra función combine_estimates el cual permite obtener los intervalos creibles (no de confianza). Pero antes, debemos realizar una pequeña transformación a nuestra anterior tabla estimacion, esto con el fin de obtener nuestros resultados desagregados por sexo (y proceder luego a combinar). Esto lo presentaremos a continuación:

tabla_sampler <- estimacion %>%
    group_by(stratum_name, replicate) %>% 
    mutate(sample_number = glue("sample_{row_number()}")) %>% 
    pivot_wider(id_cols = c("replicate", "stratum_name", "sexo", "n_obs"),
                values_from = N,
                names_from = sample_number) 

paged_table(tabla_sampler, options = list(rows.print = 10, cols.print = 5))

En primera instancia vemos que grupamos por nuestra columna stratum_name y cada una de nuestras 10 réplicas. Adicionalmente generamos una columna denominada sample_number en el que {row_number()} corresponde a cada fila de N. Es decir, articulando esto con nuestra siguiente línea de código, vemos que pivot_wider convierte nuestro dataframe de un formato largo a ancho, o, en otras palabras, los valores que están en N se “trasladan” a las columnas que comienzan por sample_. Como habíamos explicado anteriormente, cada estrato (por ejemplo 1990_1993-ETNICO-HOMBRE para la réplica (R1) contiene 1000 muestras aleatorias de la distribución posterior, por lo que cada valor de este N se expanderá a cada columna de sample_ y por tanto, esas serán 1000 columnas, desde sample_1 hasta sample_1000.

Teniendo esta transformación, procederemos a agrupar por nuestra variable de interés (sexo) y réplica, seguido de la suma de estas muestras de la distribución posterior de cada estrato que forma parte de este agrupamiento.

tabla_sexo <- tabla_sampler %>% 
    group_by(sexo, replicate) %>% 
    summarize(across(c(starts_with("sample_"), "n_obs"), sum)) 

paged_table(tabla_sexo, options = list(rows.print = 10, cols.print = 5))

En otras palabras, primero agrupamos por sexo y réplica, es decir -por ejemplo- agrupar por “HOMBRE” y la réplica 1, “HOMBRE” y réplica 2, y así sucesivamente. Cuando tengamos este tipo de agrupación (como lo vemos en la tabla_sexo) tomaremos -por ejemplo- la réplica 1 (“R1”) y la categoría “HOMBRE” y sumaremos los valores para el sample_ (es como tomar: ejemplo <- tabla_sampler %>% filter(replicate == "R1", sexo == "HOMBRE") y sumar: sum(ejemplo$sample_1), dandonos lo que nos está arrojando la primera fila de la columna sample_1, es decir, 17 588).

Sumado a esto, volveremos a transformar nuestros datos de ancho a largo, es decir, volveremos a una estructura muy parecida a tabla_sampler, pero esta vez nuestro N muestra la suma de las muestras de la distribución posterior por sexo y réplica que hicimos en el paso anterior:

sexo_group <- tabla_sexo %>%
    group_by(sexo) %>%
    pivot_longer(starts_with("sample_"), 
                 names_to = "replicate_num", 
                 values_to = "N") %>%
    ungroup() %>%
    select(-replicate_num) %>%
    rename(stratum_name = sexo) 

paged_table(sexo_group, options = list(rows.print = 10, cols.print = 5))

Para finalizar, vemos que esta información está desagregada por sexo y réplica, pero, como sabemos, deseamos ver los resultados únicamente por sexo. Por tal razón agruparemos esta información por sexo (que ahora se denomina stratum_name):

final_agrupacion <- sexo_group %>%
    group_by(stratum_name)

estimates_tabla_sexo <- final_agrupacion %>%
    group_split() %>%
    map_dfr(.f = verdata::combine_estimates) %>%
    mutate(est_lo_p = round(N_025/sum(N_mean), digits = 2)) %>%
    mutate(est_lo_p = ifelse(est_lo_p < 0, 0, est_lo_p)) %>%
    mutate(est_mn_p = round(N_mean/sum(N_mean), digits = 2)) %>%
    mutate(est_hi_p = round(N_975/sum(N_mean), digits = 2)) %>%
    mutate(est_hi_p = ifelse(est_hi_p > 1, 1, est_hi_p)) %>% 
    bind_cols(group_keys(final_agrupacion))

paged_table(estimates_tabla_sexo, options = list(rows.print = 4, cols.print = 5))

Vemos pues que, después de agrupar, tenemos los resultados derivados de estimates_tabla_sexo cuyo código permite separar nuestras categorías (hombre y mujer), luego utilizamos la función de map_dfr en el que podemos aplicar combine_estimates a cada una de nuestras categorías. Además, calculamos las proporciones para cada columna en el que nos aseguramos que las proporciones no sean menores a cero (respecto al límite inferior) y mayores a 1 (respecto al límite superior). Por último -con bind_cols- unimos estos resultados con las variables de agrupación, es decir, con sexo que la denominamos stratum_name.

Así, por ejemplo, podemos observar que el número de víctimas mujeres -menores de edad- de reclutamiento está dentro del rango de 7 326 y 9 547, con una probabilidad del 95%, siendo 8 290 el valor más probable, es decir, el promedio.